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^ ! ABSTRACT 

We investigate self-similar magnetohydrodynamic (MHD) processes in an isothermal 
| self-gravitating fluid with a quasi-spherical symmetry and extend the envelope expan- 

Qs, ■ sion with core collapse (EECC) solutions of Lou & Shen by incorporating a random 

' magnetic field. Magnetized expansion-wave collapse solutions (MEWCS) can be read- 

l/"") , ily constructed as a special case. This inside-out MHD collapse occurs at the magne- 

■ tosonic speed and magnetized EECC solutions are obtained systematically. Stagnation 

| surfaces of EECC solutions that seperate core collapse and envelope expansion propa- 

gate at constant speeds either sub-magnetosonically or super-magnetosonically. These 
similarity MHD solutions show various behaviours, such as radial inflow or contraction, 
core collapse, oscillations and outflow or wind as well as shocks. For solutions to go 
+-> ■ across the magnetosonic line smoothly, we carefully analyze topological properties of 

magnetosonic critical points. Without crossing the magnetosonic critical line, continu- 
. . , ous spectra of magnetized EECC and envelope contraction with core collapse (ECCC) 

solutions are readily obtained. Crossing the magnetosonic line twice analytically, there 
kj^j ' exists an infinite number of discrete magnetized EECC and ECCC solutions. One or 

. more sub-magnetosonic stagnation surfaces associated with these discrete solutions 

' travel outward at constant yet different speeds in a self-similar manner. In addition to 

the EECC shock solution which could change the central accretion rate, the magnetic 
field can also affect the core accretion rate. As the magnetic parameter A increases, 
the core accretion rate appropriate for the MHD EWCS becomes larger. Under the 
frozen-in approximation, magnetic fields in the envelope expansion portion would scale 
as B oc r _1 , while in the core collapse portion they would scale as B cx r -1 / 2 . We 
discuss several astrophysical applications of EECC similarity solutions to the forma- 
tion process of proto-planetary nebulae connecting the AGB phase and the planetary 
nebula, to supernova remnants, collapse of magnetized molecular cloud, to H II clouds 
surrounding massive OB stars and to a certain evolution phase of galaxy clusters. 
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1 INTRODUCTION 

The process of gravitational collapse is fundamental in var- 
ious cosmological and astrophysical contexts, and is cer- 
tainly one important phase for dynamical processes of 
star formation. Over decades, observational and theoreti- 
cal studies have been combined together to outline a plau- 
sible scenario for such gravitational collapses in gaseous 
molecular clouds, taking account of geometrical factors 
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and other main physical considerations (e.g., Bodenheimer 
& Sweigart 1968). Self-similar solutions under the simpli- 
fying assumption of spherical symmetry have long been 
pursued in the past J|Lareon|_^)69dJPejistc^^ 
[ 19771: iHunterl Il977t Il986l : IWhitworth fc Summers! Il98d 
iFoster fc Chevalierll99.C Independently, IPenstonl (|l969 albjl 
and lLarsonl (I1969allbl) started a simple analytical description 
of spherically symmetric isothermal collapse under the self- 
gravity that reproduced essential features of more detailed 
numerical simulations. In a seminal paper that introduced 
the current par adigm for a n 'inside-out collapse' in low-mass 
star formation, IShul (Il977f) derived the so-called self-similar 
expansion wa ve collapse solution (EWCS). Subsequently, 
IHunterl lll977f) carried out a more comprehensive search and 
found an infinite num ber of discrete solutions i n the 'com- 
plete solution space'. IWhitworth fc Summers! Jl985) have 
considerably expanded these solutions into a two-parameter 
continuum by allowing for weak discontinuities across the 
sonic critical line. 

For molecular clouds that are initially near a marginally 
stable equilibrium, the Larson-Penston (LP) solution pro- 
vides a sensible approximation to the condition near the 
cent r e of t he point mass (i.e., a protostar) formation epoch. 
IShul Jl977l) derived the EWCS that has a free-fall core col- 
lapse with a mass density profile of p oc r _3//2 , a radial 
infall speed profile u oc r- _1//2 , a constant central mass 



accretion rate at r — > + , and a static envelope of a 
singular isothermal sphere (SIS) with the collapsing front 
moving outward at the constant sound speed a, where r 
is the radial distance from the centre. This EWCS sce- 
nario of an in side-out collapse for st a r form ation has been 
advocated by IShu. Ad ams fe Lizanol i l987t) and compared 
with observations ([A dams. Lad^fc Shul Il987l: IZhou et alJ 



1993l:IChoi et al.ll995l:ISaito et alll999HHarvev et all200ll) . 

Foster fc Cheva lier (1993) studied the gravitational collapse 



of an isothermal sphere by hydrodynamical simulations and 
recovered the LP solution in the centr al region w here a core 
forms and the self-similar solution of IShul l|l977l) when the 
ratio of initial outer cloud radius to core radius is > 20. 

Molecular gas clouds from which many stars come into 
being are structurally complex and highly filamentary in 
many places. The origin of such complex structures is not 
well understood at present. It may be partially produced by 
regular and random magnetic fields on different scales (e.g., 
Balick & Frank 2002). The formation of filaments and sheets 
may be necessary to the fragmentation of molecular clouds 
into clumps of order of the final stellar masses, although 
the relation between the initial cloud mass and the mass of 
the resulting stars remains unclear (e.g., Hartmann 1998). 
Due to the for mation of filamentary structures in numeri- 
cal simulations JPorter et al]ll994tlKlessen fc Burkertl2000t 
lOstriker et alJl200ll) and the appearance of numerous fila- 
ments in observations (e.g., Falgarone et al. 2001), analyses 
on processes of filament formati on and evolution have also 
been pursued in the literatu re ( Kawachi^k Hanawa 1998; 
lHennebelldl2003l:lTillev fc Pudritdl2003l: IShadmehrill2005F) . 

In view of the ubiquitous presence of magnetic fields 
in the interstellar medium and their nontrivial dynamic 
and diagnostic roles in astrophysical plasma systems, dy- 
namical collapses of magnetized clouds have generated con- 
siderable interests. In some cases, magnetic fields may 
not be dynamically important, yet can offer valuable 
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diagnostic information thr o ugh c yclotron or synchrotron 
omissions. IChiueh fc Choi (1994) carried out an inves- 
tigation of dynamical collapse in nonrotating magne- 
tized gas clumps of grossly spherical shape. They fo- 
cussed on the magnetic pressure effect with the mag- 
netic tension effect being ignored. Effects of magnetic 
fields ha ve also been extensively investigated in various 
contexts ([C ontopoulos et al. 1998; Krasnopolsky^^ttn^ 
| 2002t lHennebelldl2003l : iTillev fc Pudritd 120031: IShadmehr 
2005). In terms of theoretical formulation and astrophysical 
applications, we note that magnetic field can make a tenu- 
ous collisionless gas to behave more fluid like as evidenced 
by fusion plasma experiments and solar wind observations. 
As a result, we could study dilute magnetized plasmas on 
large scales using the MHD description. 

Thermal ra dio emission s from planet ary nebulae have 
been observed IIKwokl 11982. 1985J, I1993T) and might per- 
haps indicate the presence of magnetic fields. By exten- 
sive observations of various morphologies, it is most likely 
that magnetic field should influence the shapes and shap- 
ing of PNe (e.g. Balick & Frank 2002). Although not a 
PN, the remnant morphology of a triple ring system of 
the supernova SN 1987A may also hint at important roles 
of ma gnetic field as modelled through numerical simula- 
tions llTanaka fc Washimill2002T) . Within the collapsed core 
of a PN, magnetic dynamo processes may also operate for 
some time to generate and sustain a strong magnetic field 
(e.g., ~ 10 4 — 10 9 G) associated with a proto- white dwarf 
l|Blackman et alj|200ll) . 

The role of magnetic field in our model analysis is char- 
acterized by one important parameter A as defined by equa- 
tion 11611 . Different systems of astrophysical objects involve 
various values of A. The estimated A parameter for the Crab 
Nebula is of the order of a magnitude ~ 10 5 . For a galaxy 
cluster, the estimated A is approximately 0.2. For star form- 
ing regions, the estimated A is typically 0.01. For planetary 
nebulae, the typical A is approximately ~ 0.003. 

In addition to collapses and flows with an isother- 
mal equation of state, gravitational collapses with 
polytropic or logotropic equations of state have also 
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2003; Shadmchri 2005). With a polytropic equation of state, 
the gravitational bound system have different tempera- 
tures at different radii. It is also possible to include effects 
of radiative losses in the energy equation (e.g., Boily & 
Lynden-Bell 1995; Murakami, Nishihara & Hanawa 2004). 
Stability of self-similar gravit ational collapses was also 
consi dered by several author s l|Ori fc Piranlll988l:l~ 



Hanawa 



[ 1997 1 : lHanawa fc Matsumotd Il999t l200d: lLai fc Goldreich 
hoOOULailfeOOOMLai fc Goldreich! J2000l) studied the growth 
of nonspherical perturbations in the collapse of a self- 
gravitating spherical gas cloud. They found that nonspher- 
ical perturbations damp in the subsonic region and grow 
in the supersonic region. They mentioned pote ntial ap plica- 
tions to core collapse of s u pernova exp losions llChendll978l : 
iGoldreich fc Weber! Il980l: lYahill Il983l) . where asymmetric 
density perturbations may lead to asymmetric shock propa- 
gation and breakout, giving rise to asymmetry in the explo- 
sion and a kick to a new-borne neutron star. Global stability 
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analysi s wit h a pol ytropic equation of state was also carried 
out by lLail (2000). Slow rotation effects on gravitational 
collap ses has been considered by iTerebev. Shu fc Cassenl 
(1984). As an important extension of Shu (1977), the gen- 
eral relativistic analysis of collapse or accretion onto a black 
hole has been recently cond u cted b y Cai & Shu (2005). 

Recently, iLou fc Shenl l|2004h examined self-similar 
spherical isothermal collapses and derived new solutions in 
the semi-complete solution space; this is i n reference to the 
'complete solution space ' of lHunteJ l)l977h . The gas flow at 
large r, together with the two distinct asymptotic solution 
behaviours near the origin, lead to an infinite number of 
discrete self-similar solutions that are analytically smooth 
across the sonic critical line. Those solutions characterized 
by envelope expansion with core collapse (EECC) are of par- 
ticular interest. In various contexts of astrophysical systems, 
it would be natural and important to include effects of mag- 
netic fields in a collapsing cloud. Our main motivation here 
is to study this generalized class of EECC solution in the 
presence of a random magnetic field. Meanwhile, we broaden 
and enrich all the known similarity solution structures in the 
most general manner. 

In the presence of magnetic field, the exact spherical 
symmetry for a magnetized flow no longer exists in the rig- 
orous sense. Nevertheless, morphological observations of su- 
pernova remnants, such as the Crab Nebula, show that gross 
quasi-elliptical structures do exist on large scales. With these 
empirical considerations in mind, we develop an MHD for- 
malism for a quasi-spherical symmetry for large-scale mag- 
netized flows. The close analogy between the mass conser- 
vation equation and the magnetic induction equation gives 
rise to a special relation connecting the gas mass density, 
the magnetic field component transverse to the radial flow 
direction and the radius r. With this relation, magnetic 
field can be readily incorporated into the EECC framework 
i. This paper is to investigate the mag- 
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netic field effects on the EECC solutions under the approx- 
imation of a quasi-spherical symmetry. 

Besides these shock-free globally smooth solutions, 
iTsai fc Hsul lll995l) constructed a self-similar shock solution 
describing a situation in which a central thermal or kinetic 
energy release initiates an outgoing shock during a proto- 
stellar collapse of a low-mass star. Such a shock solution, 
matched with a static SIS envelope, can have a central free- 
fall collapse, yet with a central m ass accretion rate lower 
than that predicted by the EWCS. ITsai fc Hsul (Il995l) also 
constructed a shock expansion solution with a finite core 
density matched with a static SIS envelope; this solution 
has been recently generalized and subs umed into a fam ily 
of 'champagne flow' shock solutions by llShu et al.ll2002l) to 
model expansions of H II regions surrounding massive OB 
stars after th e passage of an ini tial ionization front. To be 
more general. IShen fc Loul B2004I) studied EECC shock solu- 
tions allowing outflows or inflows far away from the central 
region. In terms of modelling proto-stellar systems, this flex- 
ibility can accommodate a variety of real possibilities. Most 
recently, Bian & Lou (2005) explored the parameter space 
more systematically to construct various isothermal shock 
solutions, including for example, twin shocks, and so forth. 

This paper is structured as follows. The basic MHD 
equations and the self-similar transformation are presented 
in § 2. In § 3, we provide a comprehensive analysis of the 



resulting nonlinear MHD ODEs. In § 4, we solve the non- 
linear MHD ODEs numerically with the help of necessary 
analytical analysis given in § 3. In § 5, we give the physical 
interpretations of our results. 



2 PHYSICAL ASSUMPTIONS AND 
BASIC MODEL FORMULATION 

We adopt the isothermal approximation for the gas under 
consideration. By the quasi-spherical symmetry of a mag- 
netized flow system, the basic MHD equations in spherical 
polar coordinates (r, 9, </>) are: 



dp 1 d(r 2 pu) 



dt 



dr 



o , 



dM 2 

— — = iirpr 
or 



(1) 
(2) 



for the mass conservation, where u is the bulk radial flow 
speed, p is the gas mass density and M(r, t) is the enclosed 
mass within r at time t. By the definition of M(r, t), we have 
the equivalent of equation Q, namely 



dM | u dM 
dt dr 







(3) 



The radial momentum equation is simply 
a 2 dp GM 1 d 



du du 
dt dr 



p dr 



8typ dr 



<B1> 



<Bj> 
4-Kpr 



(4) 



where —d&/dr = —GM(r,t)/r 2 with <E>(r, t) being the grav- 
itational potential and B±_ is the magnitude of the magnetic 
field transverse to the radial direction. The Poisson equa- 
tion relating the gas mass density p and the gravitational 
potential $ is automatically satisfied. Here in equation Q, 
we keep the mag netic tension force term that was ignored 
in equation (3) of lChiueh fc Choul l|l994h . 

The MHD energy conservation takes the form of 



d J pu 2 
dt — +ap 



In ^ -1 



8tyG [ dr ) + 8tt 



+ - 



1 d 



2 

r up 



i 2, / P 
+ a In 



1 d . 2 Id 



f>c 

r 2 $ d 2 $ 
dr V MJdrdt 



+ V 



B 

4tt 



x (u x B) 



, 



(•») 



which generalizes equation (4) of Lou & Shen (2004) by in- 
cluding the magnetic field B; the magnetic energy density 
and the Poynting flux density can be readily identified, re- 
spectively. We derive this MHD energy conservation equa- 
tion when all magnetic field components are included. 

To reach the above MHD energy conservation equation, 
we have taken the electrical conductivity to be infinite, such 
that the magnetic induction equation takes the form of 

^Vx(uxB). (6) 

The three component forms of this magnetic induction equa- 
tion arc 
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dB r 

dt 



OBg 

dt 



dt 







r sin 9 89 



sm9(u r Bg — UgB 
1 d 



r sin 9 d(j) 



(Ut/,B r — UrBtf, 



-{u g B. 



4> ' 



ld_ 

r dr 



r(u r Bg — ugB r 



(7) 



(8) 



ld_ 

r dr 



r(u$B r — u r B^j 



1 d 

- — (ugBj, - U^Bg) 

r off 



(9) 



where u r = u in our formalism. We just consider radial bulk 
flow u r 7^ by setting ug = u$ — 0, equation 10 then 
becomes 



dB r 
dt 



ii r (Vi ■ Bj 



(Bj_ • V)« r 



(10) 



Using the divergence-free condition of magnetic field V ■ B 
0, equation (1101 becomes 



1 D 

2Dt 



[(r 2 B r ) 2 ] = r 4 B r (B ± ■ V)u r 



(11) 



where D/Dt = d/dt + u r d/dr. The 9— and (^—components 
of the magnetic induction equations © and under the 
simplifying approximation of ug = and = 0, become 



D /B 6 
— - m 

Dt 



19,2 



i a 2 



(12) 



(13) 



The continuity equation of mass conservation takes a similar 
form in parallel, namely 

D 



^1»P = -^V). (14) 
From the above three equations (11211 — 1141 . we obtain 

— - = const , — - = const , — — = const , (15) 
pr pr pr 

where B± = (£>f + B 2 ,) 1 ^ 2 is the magnitude of the magnetic 
field component transverse to the radial component. 

Based on the last of integrals 1151 . we introduce a simple 
relation 



B\ = 167r 2 AGpV 



(16) 



with the magnetic parameter A being a p roportional con- 
stant. In comparison with the work of IChiueh fc Choul 
(1994), it is apparent that 



A 



f3(x) 



where 0(x) is a dim ensionless function introduced by 
IChiueh fc Choi (Il994h and a is our reduced density. Note 
that the ratio of the Alfven speed va to the isothermal sound 
speed a is given by 



VA 
a 



1/2 



and 



VA 



B ± 

(47TP) 1 / 2 



Similarity solutions with x ranging from x — * +oo to x — > + 
are referred to as the semi-complete solutions (Shu 1977; 
Whitworth & Summers 1985; Lou & Shen 2004; Shen & Lou 
2004). Similarity solutions with x ranging fro m x — » — oo to 
x — » +00 are referred to as complete sol utions dHunterll977l . 
1986). The perspe ctives of complete (iHjrrrkej ll9Jjl Ll2£fl) 
and semi-complete <ShJl977tl\Vhitworth fc Summersl^g!! 
lLou fc Shenll2004ft similarity solutions are both valid with 
appropriate physical interpretations and with the proper 
specification of an initial moment. For the invariance under 
the time reversal transformation of self-gravitational mag- 
netofluid equations Q — QJ, namely, t — > —t, p — > p, 
u — > —u, the semi-complete and complete similarity so- 
lutions are in fact closely related to each other. Complete 
Hunter type solutions can be decomposed into two branches 
of similarity solutions in the semi-complete space (Lou & 
Shen 2004). Futhermore, the time reversal transformation 
can also lead to different interpretations and application of 
self-similar solutions. In this paper, we focus on the semi- 
complete solutions. For ideal MHD equations with an infi- 
nite electric conductivity, the magnetic induction equation 
© is invariant under t — ► —t and v — > — v transformation. 
In fact, the corresponding transformation of either B — > B 
or B — *• — B is allowed. If we take the magnetic diffusion 
term into account, these invariant properties are no longer 
valid. 

We now introduce the self-similar transformation: 



p(r,t) = 



a(x) 
47rGt 2 



M(r,t) = -g-m(x) 
u(r, t) = av(x) , 
*(r,t) = a 2 (f>(x) , 



(17) 



S, = 



-b(x) 



VGt 

and immediately get 
b(x) = vXax , 

where x = r/(at) is the independent self-similar variable; 
the dimensionless a(x), m(x), v(x), <j>(x) and b(x) are the 
reduced dependent variables for the gas mass density, the 
enclosed mass, the radial flow speed, the gravitational po- 
tential and the transverse magnetic field, respectively. They 
ar e all functions of x onl y. Note that b{x) is related to @(x) 
of lChiueh fc Choul dl994l) by b 2 = f3. After a direct substitu- 
tion of the above transformations into the nonlinear MHD 
partial differential equations, we obtain the following non- 
linear ordinary differential equations (ODEs). 

It follows that equations @ and ^ reduce to 



m = x a 



(v — x)m + m — , 



(18) 
(19) 



where the prime denotes the derivative with respect to x. 
Combining these two equations, we immediately obtain the 
following two equations 

m = (x — v)x 2 a 
d 



dx 



(20) 
(21) 
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By equation I20H . the physical condition of a positive mass 
m(x) > is transformed into the condition of x — v > 0. 
Therefore, solutions of v(x) must lie to the upper-right of 
the straight line x — v = in the plane —v(x) versus x. 

The continuity equation of mass conservation becomes 



v + (v — x)- 



'-(x — v) 



and the radial momentum equation becomes 



2 Q 

(v — x)v + (1 + Xax ) — = —a(x — v) — 2Xxa . 
a 

Rearranging these two equations above, we arrive at 



(22) 



(23) 



(x — v) — (1 + Xax ) 



(x — v) 



a(x — v) 



(24) 



and 



(x — v) 2 — (1 + Xax 2 ) 



a 
a 



(x — v) 



X 



+ 2Xxa . (25) 



Note that by setting A = for the absence of magnetic field, 
we readily recover the formulation of Lou & Shen (2004). 
For smooth and regular solutions crossing the magnetosonic 
line, the critical conditions are 



(x — v 
(x 

(x — V 



a(x ■ 



12 



= , 



/')" - (1 + Xax ) = 
2 



a 



X 



'-(x — v) 



+ 2Xxa = 



(26) 
(27) 
(28) 



in order to determine finite first derivatives v' and a', re- 
spectively. In fact, only two of the above three requirements 
are independent. For i-jj / 0, the above three equations 
are equivalent to the following three equations 

2 



(29) 

(30) 
(31) 



2 2 . , 3 4 . 

a x + Xa x — 4 , 

(x - v) + 2Xx = (x- v f . 

For A = 0, equation 1311 becomes (x — v) = ±1 and equation 
<■'-)! H reduces to a = 2/x. Given the positiveness of x — v > 0, 
we then come to t he familiar critical condition x — v — 1 
jLou fc Shenl2004l) . Equation gfl for (x-vf -(l + Xax 2 ) = 
is equivalent to the condition (a; — v) 2 — (a 2 + v\)/a 2 
in dimensional form, directly related to the magnetosonic 
condition. 



3 ANALYSIS OF NONLINEAR ODES 
3.1 Special Analytical Solutions 

With (f> = ma; -2 , an exact global analytical solution of non- 
linear ODEs O and appears in the form of 



2 d(j> 2 2X 1/2 , . 

w = 0, a = — , m = 2x, = -, b= , (32) 

x A ax x x 

which describes a magnetostatic equilibrium for a quasi- 
spherical system lEberdll955t lBonnorlll955 IChandrasekharl 
Il957t IShulll9T7fl . In physical dimensions, we simply have 



P = 



2irGr 2 



M = 



2a 2 r 
IT 



= ,Xa* 



Gr 2 



In this solution, the magnetic pressure gradient force den- 
sity and the magnetic tension force density cancel each other 
(i.e., a nearly force-free magnetic field), that is, their com- 
bined effects on the entire system are absent. To be more 
specific but in dimensionless form, the sum of the magnetic 
pressure gradient and tension force densities is simply 

-X(a'x 2 + 2ax) = . 

Another exact global analytical similarity solution of non- 
linear ODEs lO and J23 is 



2 

v=-x 



2x 3 



(3 + 18A) 



3(3 + 18A) ' 



2x 



dx 3(3 + 18A) 



2X 1 ' 2 x 
(3 + 18A) 



(33) 



In physical dimensions, we simply have 
1 » , 2r 3 



P - 



2(3 + 18A)ttG£ 2 ' 

_ 4Ar 2 

~ G(3 + 18A) 2 t 4 ' 



M = 



3(3+ l8X)Gt 2 



Within a finite radius r, the transverse magnetic field 
strength \B±\ scales as t~ 2 . This solution passes the mag- 
netosonic critical line at x = 3(1 + 6A) 1,/2 as would be evi- 
dent in our later discussion about the MHD generalization 
of the Larson-Penston-type solutions. In the absence of mag- 
netic field wit h A = 0, this solution passes t he sonic critical 
line at x = 3 dWhitworth fc Summerslll985h . Such solution 
corresponds to a nonrelativistic Einstein-de Sitter expan- 
sion with magnetic field. Here, the 'cosmic mass density' is 
p = l/[(6 + 36A)-7rG£ 2 ] and the relevant 'Hubble constant' is 
H = 2/[3t(l + 6A) 1/2 ] that decreases with time t. 

Another exact global solution is the singular solution 



2 

ax 



and [x — v) 2 — (1 + Aaa; 2 ) = 



(34) 



which actually specifies the magnetosonic critical line. This 
set of equations for the magnetosonic critical line shall be 
discussed in more details presently. 

The second condition of equation l|34^ is equivalent to 

a(x — v) = (v 2 A + a 2 ) 1/2 ; 

the right-hand side is the magnetosonic speed and the left- 
hand side is the wave front speed relative the local flow 
speed. When the Doppler-shifted flow speed coincides with 
the local magnetosonic speed, the singularity would arise. 
The situation appears to be similar to the discussion of the 
'nozzle condition' in the steady solar wind (IParkerlll963h . or 
to the discussion o f the critical accretion rate in a steady 
spherical accretion tf3ondilll952l) . 
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3.2 Results for Crossing the Magnetosonic 
Critical Line 

The coupled nonlinear MHD equations are singular along 
the magnetosonic critical line. Because of this, we cannot 
directly obtain similarity solutions that go across the magne- 
tosonic critical line smoothly by a simple-minded numerical 
integration backward from x — > +00. To obtain solutions 
satisfying magnetosonic critical conditions 1291 — 1311 . we 
need to first derive non-trivial first derivatives of v and a by 
applying the L'Hospital rule along the magnetosonic critical 
line. Higher-order derivatives of v and a across the mag- 
netosonic crit i cal line can also be determin ed accordingly 
jHunterlll98r3 : IWhitworth fc Summersl lT985^ . It is straight- 
forward to derive a quadratic equation for z = v , namely 



(v — x) 



(v 



z + 2(x — v)z - 



2 4 = o 



(35) 



(see Appendix B for a detailed derivation). 

Solving this quadratic equation, together with equation 
12911 . 13111 and 12211 . we obtain the first derivative v' along 
the magnetosonic critical line. Once the first derivative v' 
is known, we immediately obtain the first derivative a' us- 
ing equations 1221 or 12311 . By this procedure, the problem 
of crossing the magnetosonic critical line smoothly and an- 
alytically has been completely solved by means of an ef- 
fective numerical scheme. For quadratic equation 1351 . an- 
alytical expansions could be found once the coefficients of 
this equation can be given explicitly. We have solved the 
magnetosonic critical line equation as described later in this 
section, that is, we know the the value of v for a given value 
of x along the magnetosonic critical line. Thus analytical 
eigensolutions of equation 1351 can be given, as all the coef- 
ficients of the equation 1351 can be expressed in terms of x. 
However, such explicit eigensolutions of equation 1351 would 
be very complex and it would be more convenient to derive 
the two eigensolutions numerically. 

In compar i son with equations (16) — (18) of 
IChiueh fc Choi lll994T) . we find that their results about 
the magnetosonic critical point is a very special case of 
ours. By setting v = , we get similar results as those of 
IChiueh fc Choul (Il994h . Their equations (16) - (18) are 
results about the zero point x c = (1 + 2A) 1 / 2 . Expressing 
our results related to the zero point in terms of their similar 
notations, we have 



A 



Vc = , v cl 



2(xl - 1) 



X c 

2x1 
3a; 2 - 1 



, «c2 = , 



4(2^ - 1) / _4_ 
Qcl x%(3xc — 1) ' &c2 x% ' 

n o\ ' 2 . ox 2 8A(:r 2 - 1) 



x 3 c{3x 2 c - 1) ' 



(3 c2 = 2Aa c Q c2 Xc + 2\a 2 c x c = — ^-rr 

xt 



(36) 
(37) 
(38) 
(39) 



Note that A > and thus x c ^ 1 iChiueh fc Choulll994) . 
The derivative of b = X 1 ' 2 ax can be calculated once re- 
sults of a are known. The key difference lies in the fact that 
our momentum equation contains the magnetic tension force 



term that was ignored in their formulation. These results are 
crucial for constru cting MHD expansion-wave colla pse solu- 
tions (MEWCSs) (IShulll977l:IChiueh fc CI10U1I19941) . 1 

By equation 1351 . there are two types of eigensolutions 
for z = v' at each point along the magnetosonic critical 
line. When the two roots of equation 1351 are of opposite 
signs, type 1 and type 2 eigensolutions are referred to those 
with negative and positive roots of the equation 1351 . re- 
spectively. When the two roots of equation 1351 are of the 
same sign, type 1 and 2 eigensolutions are those with smaller 
and larger absolute values respectively, following the nomen- 
clature of Shu (1977) and Lou & Shen (2004). In Lou & 
Shen (2004), type 1 corresponds to dv/dx = 1 — 1/x, and 
type 2 corresponds to dv/dx = 1/x*. In the open range of 

< x < 2, type 1 and type 2 are exactly defined as such. 
When x > 2 (relevant to the LP-type solution), dv/dx of 
type 1 has a larger absolute value, while dv/dx of type 2 has 
a smaller absolute value, that is, their magnitudes reverse 
for nodal points. In our current definition, no magnitude re- 
versal would happen. When the point is a nodal point, type 

1 and type 2 solutions remain always the smaller and larger 
ones respectively for the absolute value of dv/dx. To summa- 
rize, our definition is not defined by the explicit expressions 
of dv/dx (e.g., 1 — 1/x, and 1/x*), but by their magnitudes 
and signs. It is easier to keep in mind their relevant physical 
properties. 

As a necessary consistent check, we note that when A = 
and x — v — 1, the familiar results of Shu (1977) are readily 
recovered, namely 



2z + 2z - 



2(a:-l) 



= 



(40) 



The two eigensolutions of this quadratic equation are simply 
1 . 1 



Z2 = l 



[see the Appendix of Shu (1977)]. 

Expanding equation 1)31^ and regrouping terms, we ob- 
tain the following cubic equation in terms of v{x) 
3 



v + a,2V + aiv + oo = 



(41) 



describing the magnetosonic critical line, where the three 



coefficients are 0,2 



-3x, a\ — 3x — 1 and ao = 2Xx + x - 



x . This equation can be solved for v in terms of x. Let zi, 
Z2, Z3 be the three roots of cubic equation 141fl in terms of 
x. By the properties of a cubic equation, it follows that 



21 + Z 2 + Z 3 — —0,2 , 
212:2 + 22Z3 + 2321 = 0,1 

212223 = — ao ■ 



(42) 
(43) 
(44) 



Following the standard procedure of solving a cubic alge- 
braic equation, we then introduce two dimensionless param- 
eter q and r such that 



ai a 2 

(aia 2 - 3a ) a% 



6 



27 



(45) 
(46) 



1 We note in passing that equations (1) — (3) and (8) — (10) of Chi- 
ueh & Chou (1994) are somehow inconsistent. The magnetosonic 
critical condition in their analysis is correct. 
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It then follows that 



3 , 2 

q + r 



(4a i - a 2 ) 



+ ao 



+ 



aoaia2 



(47) 



For g 3 +r 2 > 0, the cubic equation possesses one real root for 
v and a pair of complex conjugate roots for v. For q 3 +r 2 — 0, 
all three roots of v are real and at least two of them are 
equal. For g 3 + r 2 < 0, all three roots of v are real. In cubic 
equation 1411 . v is the unknown for a given x. In solving for 
v, we express the three roots z% as 



zi = (si + s 2 ) ~ a 2 /3 
(si + s 2 



22,3 = 



fl2 . iV3, 
-J ± — (S1 



«2 



(48) 
(49) 



where the two parameters sj and s 2 are defined by 

_ r i/3 i 2n1/2i1/3 

Si, 2 = [r ± (q + r ) J . 

Specific to the current problem, we simply have 
q =-1/3, (50) 
r = -xX (51) 
and 

g 3 +r 2 = a; 2 A 2 - 1/27 , (52) 
when 

x > 1/(3V3A) ; (53) 

the only real root is the one given by equation (1481 . Given an 
explicit expression of this root, we can analyze the nodal or 
saddle properties along the magnetosonic critical line (e.g., 
Jordan & Smith 1977). For x < (3V3X)" 1 , we need to iden- 
tify the suitable root lying in the required parameter regime 
to satisfy x — v > 0. This is somewhat more involved than 
the case of x > (3V3A) -1 . 

For the case of x ^ (3\/3A) _1 , we adopt an alternative 
solution procedure by introducing a transitional variable 



(3y/3xX) 



(54) 



where 8 lies in the range 0^8^. ir/2. The three roots of the 
cubic equation can be express ed more compactl y and neatly 
in terms of this 9 parameter iPress et al.lll986T) . namely 



z\ 



22,3 



-q cos 



2^3 



-2y/—q cos 



2^3 



\ C12 






J 3 






e\ x 

3j + 3 




(55) 








e± 2tt^ 


3 




e± 2tv\ 

3 ) 


X 

+ 3' 


(56) 



Since 8/3 falls in the range of ^ 8 ^ vr/6, it follows that 
cos(0/3) is positive while cos(S/3+27r/3) and cos(0/3— 27r/3) 
are both negative. Therefore, Z2 and 23 must be positive, 



Magnetosonic Critical Line 
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^=1.0 







-3 -2 -1 



Figure 1. The magnetosonic critical lines for A = 0.1, 1/6, 1.0 in 
upper, middle, lower panels, respectively; the vertical coordinate 
is the negative reduced radial speed —v = —u/a and the horizon- 
tal coordinate is the independent similarity variable x = r/(at). 
In the upper panel for A = 0.1, the zero point is to the left of the 
turning point. In the middle panel for A = 1/6, the zero point 
and the turning are at the same location. In the lower panel for 
A = 1, the zero point is to the right of the turning point. This 
diagram is consistent with our analytical analysis. The diagonal 
straight line from the upper-left to the bottom-right in each panel 
is x — v = 0. 



while 21 may be either positive or negative. The zero point 
(i.e., —v — 0) of the magnetosonic critical line is at 

x = (1 + 2A) 1/2 . 



The other special location is 
x = 1/(3V3X) 



(57) 



(58) 



that separates — v regimes for one root (for larger |x|) and 
three real roots (for smaller |a;|). By equating these two spe- 
cial locations, we have 



A, = 1/6 . 



(59) 



We refer to the point with abscissa x = l/(3y3A) in the 
fourth quadrant as the turning point along the magnetosonic 
critical line. In Figure 1, we see that for A = 1/6, the zero 
point and the turning point of the curve in the fourth quad- 
rant have the same abscissa x as expected. For A < 1/6, 
the abscissa of the zero point is smaller than that of the 
turning point in the fourth quadrant, while for A > 1/6, 
the abscissa of the zero point is larger than that of the 
turning point of the curve in the fourth quadrant. From 
the three panels in Fig. 1, we conclude that the magne- 
tosonic critical lines are connected smoothly by the two so- 
lutions 1551 and 1481 at x = (3\/3A) _1 . A quick summary for 
the magnetosonic critical line in terms of —v is as follows: 
when x ^ (3\/3A) -1 , z\ = — (2/v / 3) cos(#/3) + x/3; when 
x ^ (3\/3A) _1 , z\ — S\ + S2 + x/3. Once the magnetosonic 
critical line in terms of v(x) is obtained, the corresponding 
magnetosonic critical line for a(x) can be readily determined 
by equation 1291 . In Figure 2, we explore variations of the 
magnetosonic critical lines for different values of A. By the 
time reversal symmetry of the magnetosonic critical line, we 
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Variation of Magnetosonic Critical Line with X 




Figure 2. We choose A = 0.1, 0.2, 0.3, 0.5, 1.0, 2.0, 5.0, 10.0 
(from the lower part to the top part in order in the first and 
fourth quadrants); the vertical coordinate is the negative reduced 
radial speed —v = —u/a and the horizontal coordinate is the 
independent self-similarity variable x = r/(at). It is clear in this 
figure that in the first and fourth quadrants, all the magnetosonic 
critical lines lie above the straight line —v = \/3/3 — x. This 
information is essential for our analysis of the nodal and saddle 
points in Appendix B. The diagonal straight line from the upper- 
left to the bottom-right is x — v = 0. 



plot them in the entire plane. Our focus is mainly on the 
first and fourth quadrants. 



3.3 Asymptotic Solutions of Nonlinear ODEs 

Solution behaviours of nonlinear ODEs 1241 and 1251 at in- 
finity (x — > +oo) show effects of magnetic field (i.e., A / 0) 
in the higher order terms, as seen clearly from the following 
expressions of asymptotic expansions 

T7 2-A V 
v = V + h -it 

x ar 



+ 



(A/6 - l)(A - 2) + 2V 2 /3 + A{2 - A)\/3 



x" 



(60) 



A_ A(2 - A) A(4-A)V 
a ~ x 2 + 2^4 + 3^ + 



A(A - 3) (4/2 — I) — (A — &)AV 2 /4 + (2 - A)A 2 \/4 



+ 

(3 — \ax 
+ 



2J _ A^X A l (2 — A)\ A 2 (2- Af\ 
x 2 x 4 4x 6 



(61) 



(62) 



where V, A and A are three integration constants referred 
to as the speed, density and magnetic field parameters, re- 
spectively. The magnetic field parameter A comes into play 
fairly late in the expansion of large x. In the reduced velocity 
expansion, magnetic fields first enter in the x~ 3 term. In the 
reduced density expansion, magnetic fields first enter in the 
x~ 6 term. So the MHD solution behaviours at large x are 
only slightly modified by the presence of magnetic field. The 



magnitude of transverse magnetic field B± at large distances 
scales as 

„ 2 Xa 4 A 2 



Gr 2 



We may get back to nonlinear ODEs 1241 and 1251 to quali- 
tatively see the leading expansion behaviours of the solutions 
16U1 and 1611 . In the limit of x — > +oo, the leading term of 
v is the constant speed parameter V and the leading term 
of a is Ax~ 2 . By substituting the leading terms into the fol- 
lowing expressions for various relevant forces, we obtain the 
pressure gradient force in the limit of x — > +oo as 

— a /a = 2/x , 

the gravitational force in the limit of x — > +oo as 

< s A 
— a(x — v) = , 

x 

and the magnetic pressure and tension forces together in the 
limit of x — > +oo as 

~\(a'x 2 + 2ax) = . 

From these three expressions, we know that the sum of mag- 
netic pressure and tension forces actually does not affect the 
dynamic behaviour at very large x (i.e., a nearly force-free 
condition). When A > 2, the inward gravitational force over- 
powers the pressure gradient force. When A = 2, the grav- 
itational force is barely balanced by the pressure gradient 
fo rce. The lat ter situation is closely pertinent to the EWCS 
of lShl (Il977l> . 

In the limit of x — > + , the asymptotic solutions for the 
central free-fall collapse is 

v -> -2F/x 1/2 - -^x 1/2 lax- 2Lx 1/2 + ■■■ , (63) 



4F 



F/x 3/2 - ± x- 1/2 lnx-Lx~ 1/2 + 



8F 

a \ 2 2 ^F 2 

p = \a x = 

x 



+ ■ 



(64) 



(65) 



where F > and L are two constant parameters. We may 
take F = mo/\/2 with parameter mo being a measure for 
the central mass accretion rate (Lou & Shen 2004). In di- 
mensional form, the transverse magnetic energy density is 



Bl 
8n 



is 



16vrGtr ' 

The pressure gradient force in the limit of x — -> 4 

_a_ _ 3_ 

a 2x ' 

the gravitational force in the limit of x — > + is 
2F 2 

-a[x-v) = —^- , 

and the sum of the magnetic pressure and tension forces 
together in the limit of x — > + is 

' 2 AF 
-A(q x +2ax) = -1^2 ■ 

In the limit of a; — > + , the most dominant force is the 
gravitational force leading to a gravitational core collapse. 
The pressure gradient force is the second dominant one. In 
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comparison, the sum of the magnetic pressure and tension 
forces is inward and has the weakest effect as x — * + . 

The other asymptotic solution for a core expansion as 
x — > + is given by 



2 (2-ao-iaDA) x a 

3 135 

„ D(2 -3D- 18DX) 2 
D H -x + ■ 



f3 — Xa 2 x 2 



18 

XD 2 x 2 + 



(66) 

(67) 
(68) 



which is a fairly straightforward MHD generalization of the 
Larson-Penston type solutions. 

The pressure gradient force in the limit of x — » + is 

a (2 - 3D - 18DX) 

= „ 2 , 

a 9 

the gravitational force in the limit of x — > + is 
-a(x-v) = -—x , 

and the sum of the magnetic pressure and tension forces in 
the limit of x — > + is 

-\{dx 2 + 2ax) = —2XDx . 

By comparing the three forces, it is clear that the magnetic 
field plays an important dynamical role in the core expan- 
sion. This involves a radially inward net magnetic force. By 
setting A equal to zero, we readily recover the results of 
IShu et alJ J2002D . From equation we know that v > 0, 
v > and (v— x) < 0. The total time derivative of u = av(x) 
is nevertheless negative and and the magnetic energy density 
is simply 



Bi 
8tt 



XD r 



87rGt 4 

Please note that the straight line x — v = is another spe- 
cial demarcation line of importance, especially for the phys- 
ical requirement of a positive m(x). In the presence of mag- 
netic field, the magnetosonic critical lines become curved in 
a characteristic manner. In contrast, the sonic critical line 
in the absence of magnetic field is a straight line x — v — 1 
in the first and fourth quadrants. When the magnetic field 
is included, the magnetosonic critical lines become curves in 
the first and fourth quadrants. In order to go across the mag- 
netosonic critical line smoothly, the primary and secondary 
directions of the eigensolutions need to be reanalyzed ac- 
cordingly. The topological characteristics of the saddle and 
nodal points have also been changed by the inclusion of mag- 
netic field (see Appendix A for details) . 



4 NUMERICAL SOLUTIONS OF 
THE NONLINEAR MHD ODES 

We derive global solutions in the semi-complete solution 
space from x — > +oo to x — *■ + . In fact, there exists a 
one-to-one correspondence between the semi-complete and 
complete solutions, because of the invariance property for 
the time reversal symmetry discussed earlier (see figure 6 
of Lou & Shen 2004 for examples). We shall solve below 
the coupled nonlinear MHD ODEs for similarity solutions 
making use of the magnetosonic critical line analyses and 



the asymptotic solution properties presented in § 3. Numer- 
ically, we can construct various types of solutions that may 
or may not cross the magnetosonic critical line. Integrating 
from x — > +oo with three constant parameters A, V and A 
specified in asymptotic solutions 1601 and 161L we readily 
obtain those similarity solutions without crossing the mag- 
netosonic critical line. 

For global smooth self-similar solutions satisfying crit- 
ical conditions (1291 — ll-i 1 1 along the magnetosonic line, we 
need to determine the first derivatives of v and a using the 
L'Hospital rule. We then integrate leftward (backward) and 
rightward (forward) from the magnetosonic critical point. 
Based on the root properties of quadratic equation 1351 
for eigenvalues of the radial speed gradient v' , we can fig- 
ure out the nature of the magnetosonic critical points (e.g., 
Jordan & Smith 1977; Whitworth & Summers 1985). The 
magnetosonic critical line with x — v > covers the first 
and fourth quadrants. Points along the magnetosonic criti- 
cal line in the first quadrant [0 < x < (1 + 2A) 1//2 ] are saddle 
points, while points along the magnetosonic critical line in 
the fourth quadrant [(1 + 2A) 1//2 < x < +oo] are nodal points 
(see Appendix A for details). In the first quadrant, type 1 
eigensolutions correspond to negative v gradients and type 
2 eigensolutions correspond to positive first derivatives of 
v. In the fourth quadrant, type 1 eigensolutions correspond 
to lower absolute values of first derivatives of v and type 2 
eigensolutions have larger absolute values of first derivatives 
of v with both v' being positive (see Appendix A) . Numeri- 
cal integrations away from a saddle point tend to be stable. 
Unstable or neutrally stable integrat ions may occur when in- 
tegrating away from a nodal point llWhitworth fc Summers! 
119851: lHunt"erj|l986l : [Jordan fc Smitblll977l) . In the following 
we explore numerically various types of solutions that may 
or may not intersect the magnetosonic critical line. 



Global solutions without crossing the 
Magnetosonic Critical Line 



4.1 



The MHD generalization of EWCS of|Shu| (119771) can be ob- 
tained by joining two solution segments. One portion for a 
central MHD free-fall collapse is to the left of x = (1 + 2A) 1/2 
and the other for a magnetostatic quasi-spherical envelope 
is to the right of i = (1 + 2A) 1//2 . The inner part of 
MHD EWCS can be obtained by integrating leftward from 
x = (1 + 2A) 1 / 2 with the first derivatives of v and a given 
by v c i and a c \ in equations 1361 and 13 7H . respectively. Note 
that first derivatives v c2 and a c2 cannot be used to integrate 
leftward. The right portion of the MHD EWCS is the spe- 
cial analytical solution l^. In Figure 3, the MHD EWCS 
solution is plotted in a heavy solid line. This solution can 
also be obtained by taking the limit of A — > 2 + . Note that 
the MHD EWCS is tangential to the magnetosonic critical 
line at the point x = (1 + 2A) 1 / 2 . At this point, a = 2x~ 2 
and thus 

(a 2 + v 2 A )/a 2 = 1 + Xax 2 = 1 + 2A . 

In other words, we have r/(at) = (a 2 + v 2 l ) 1 ^ 2 /a correspond- 
ing to an MHD collapse wave front travelling outward at the 
magnetosonic speed (a 2 + w^) 1 ' 2 (i.e., the fast MHD wave 
speed perpendicular to the local magnetic field). 
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Figure 3. As examples of illustration, we choose A = 0.1, V = 0, 
A = 2.0, 2.2, 2.4, 2.6, 2.8, 3.0 (from bottom to top in order 
along the solution curves); the vertical coordinate is — v and the 
horizontal coordi nate is x. Th is figure is qualitatively similar to 
figure 1 of iLou &: Shenl 120041) . With a magnetic field included, 
the critical line has a different shape from that of the pure hy- 
drodynamic case. Dashed line is x — v = 0. Dot-dashed line is the 
curved magnetosonic critical line. 



In star forming regions, e.g., the B335 cloud system, a 
typical mass density in molecular clouds is of the order of 
10 -18 g/cm 3 and a typical magnetic field strength is around 
100/iG. The Alfven speed va is approximately 0.4 km s _1 , 
which is comparable to the typical sound speed a in molec- 
ular cloud cores with temperature of a few tens of Kelvin 
degrees. Therefore the inclusion of a tangled magnetic field 
would increase the central mass accretion rate by a factor of 
~ 3 or so in the scenario of a magnetized EWCS as compared 
to the purely hydrodynamic EWCS. The overall dynamical 
timescale of the infalling envelope would be shortened to a 
few times 10 yr. 

The reduced mass density a(x) is continuous at the 
'kink' point for the magnetosonic collapse front. In the case 
of an MHD generalization of Shu (1977), the central accre- 
tion rate changes with the variation of A. As A takes higher 
values, the core accretion rate becomes greater. 

As x approaches infinity, solutions should gradually con- 
verge to the asymptotic solutions I60H and 1611 . while as 
x approaches the origin, global solutions should satisfy so- 
lutions 1631 . 164H . 166H and 16711 . Those solutions without 
crossing the magnetosonic critical line can be readily ob- 
tained by numerically integrating from the asymptotic so- 
lutions 1601 and 1611 at a sufficiently large x. For those so- 
lutions crossing the magnetosonic critical line analytically, 
the information about the magnetosonic critical line and the 
gradients of v and a at the critical point is crucial for con- 
structing the solutions. 

In Figures 3 — 5, we integrate the coupled nonlinear 
ODEs 1241 and 1251 from a large enough x with V, A and 
A specified and show global solutions without crossing the 
magnetosonic critical line. In Figure 3, we numerically in- 
tegrate the equations with V = (the so-called asymptotic 
breeze solutions; Shen & Lou 2004). For V < and V > 
corresponding to inflows and outflows at x — * +oo respec- 
tively, self-similar solutions can be readily constructed. For 



0.55 
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Figure 4. We adopt A = 0.1 as in Figure 3. The corresponding 
AlfVen-to-sound speed ratios va/u of similarity solutions with 
V = 0, A = 2.0, 2.2, 2.4, 2.6, 2.8, 3.0 (from bottom to top) are 
displayed; the vertical coordinate is the speed ratio va/o. and the 
horizontal coordinate is x. 
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Figure 5. Examples of radial inflows at large x with A = 0.1, 
V = —2, A = 0.5, 0.8, 1.1, 1.4, 1.7, 2.0, respectively (from 
bottom to top); the vertical coordinate is — v and the horizontal 
coordinate is x. The straight dashed line is x — v = 0. Dot-dashed 
line is the magnetosonic critical line. 



a specified value of speed parameter V, the mass parameter 
A should meet certain criteria in order to have global so- 
lutions. Due to the nonlinearity of these ODEs and in the 
presence of magnetic field, this problem is fairly complex and 
solut ions are explored through trial and error iLou fc Shenl 
120041) . 

We now summarize some results of our numerical ex- 
ploration. For numerical solutions with V = 0, condition 
A > 2 is necessary for self-similar solutions to exist without 
crossing the magnetosonic critical line in the first quadrant. 
For V < 0, the lower limit of A for possible similarity solu- 
tions without crossing the magnetosonic critical line is less 
than 2. In our illustrating examples with A = 0.1, the lower 
limit of A is ~ 0.39. For A < 0.39, similarity MHD solutions 
would crash onto the magnetocritical line. For V > and 
A = 0.1, the minimum value of A for similarity MHD so- 
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Variations of Accretion Rate with A 



A 





0.1 


1.0 


2.0 


mo 


0.98 


1.04 


1.42 


1.73 



Table 1. Central mass accretion rates of MHD expansion- wave collapse solution (marked by the heavy solid line in Fig. 3) matched with 
an isothermal magnetostatic envelope. 



0.45 1 , 1 , 1 1 — I 1 




Figure 6. Here A = 0.1 as in Figure 5. The correspond- 
ing speed ratios of v^/a for solutions with V = —2, A = 
0.5, 0.8, 1.1, 1.4, 1.7, 2.0 respectively (from bottom to top); 
the vertical coordinate is the speed ratio va/o, and the horizontal 
coordinate is x. 



Figure 8. For A = 0.1 as in Figure 7, we display the correspond- 
ing speed ratios v^/a for solutions with V = 2 and A = 6, 8, 10, 
respectively (from bottom to top); the vertical coordinate is the 
speed ratio v^/a and the horizontal coordinate is x. 



Magnetized EECC Solution without crossing critical line 
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Figure 7. For A = 0.1, V = 2, and A = 6, 8, 10, respectively 
(from bottom to top), we readily obtain MHD EECC similarity 
solutions without crossing the magnetosonic critical line; the ver- 
tical coordinate is — v and the horizontal coordinate is x. The 
stagnation point with — v = travels outward with a super mag- 
netosonic speed. The straight dashed line is x — v = 0. The dash- 
dotted curve is the magnetosonic critical line. 

lutions without crossing the magnetosonic critical line now 
becomes generally greater than ~ 2, and varies for different 
values of V. In our specific computations with A = 0.1 and 
V = 2, the minimum value of A is ~ 5.6. For A < 5.6, the 
MHD solutions would crash onto the magnetocritical line. 



4.2 Analytical Solutions for Crossing 
the Magnetosonic Critical Line 

4-2.1 Type 1 Similarity Solutions 

We first describe the type 1 similarity solutions. By integrat- 
ing leftward and rightward from a point along the magne- 
tosonic critical line with a negative first derivative of v(x), 
we obtain the type 1 self-similar solutions. These solutions 
of v(x) i n th e first quadrant, referred to as the plus solu- 
tions by IShul <ll977l) . approach negative constant values of 

V as x — *• +oo. Following the same procedure, we can also 
construct type 1 solutions which pass through the magne- 
tosonic critical line in the fourth quadrant. These solutions 
have the properties that v(x) approaches a positive constant 

V at large x. We shall see that a class of discrete yet infinitely 
many type 1 solutions matches finite solutions I66H and 1671 
near the origin x = 0. This type of solutions were first de- 
rived by iHunteil dl977t) . Our analysis generalizes these self- 
similar solutions by including a tangled magnetic field. For 
further generalizations of LP solutions with various shocks, 
the interested reader is referred to the recent work of Bian 
& Lou (2005). 

The plus solutions of Shu (1977) (or type 1 solutions, 
see solution curves in Figure 6) with positive d(—v)/dx 
in the first quadrant have their counterparts in the pres- 
ence of magnetic field. They pass through the magnetosonic 
critical line and approach a constant negative speed V 
at large x. These solutions may be matched with an ap- 
proximate quasi-magnetohydrostatic state. As x — ► + and 
v — *• 0, we derive the following equation governing the quasi- 
magnetohydrostatic state, namely 
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Figure 9. For A = 0.1, we show two classes of examples for 
type 1 solutions smoothly crossing the magnetosonic line in the 
first and fourth quadrants, separately. The critical points are 
x t = 0.5, 0.6, 0.7, 0.8, 0.9, 1.4, 1.6, 1.8, 2.0, respectively. The 
straight dashed line is x — v = and the dash-dotted curve is the 
magnetosonic critical line. 
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Figure 10. For A = 0.1 as in Figure 9, the corresponding Alfven 
to sound speed ratios vj± /a of type 1 solutions that smoothly cross 
the magnetosonic critical line in the first quadrant. The chosen 
critical points are x* = 0.5, 0.6, 0.7, 0.8, 0.9, 1.4, 1.6, 1.8, 2.0. 



and Lou & Shen (2004). Choosing a meeting point x m and 
from different first critical point x*(\) [x*(l) < (f + 2A) 1 / 2 
and x m > x*(l)], we integrate the coupled nonlinear MHD 
ODEs forward towards the meeting point x m using the type 
2 eigensolution. With every a;*(l), we obtain a pair {v, a} 
at the meeting point x m in the so-called phase diagram of v 
versus a. For a series of a:*(l), we obtain a series of {v, a} 
pairs and thus a curve in the phase diagram. Similarly, we 
can integrate the coupled nonlinear MHD ODEs from dif- 
ferent x,(2) [x*(2) > (1 + 2A) 1/2 and x«(2) > x m ] backward 
towards the same meeting point x m to obtain another curve 
in the {v, a} phase diagram. When we integrate from x»(2) 
towards the meeting point x m , both type 1 and type 2 so- 
lutions may be chosen to construct such solutions, while 
integration from x,(l) necessarily requires the type 2 solu- 
tion. Regarding the stability of numerical integrations, we 
note that to go from a;* (2) towards the meeting point x m 
using type 1 solutions, the integration is away from a node 
along the secondary direction and is neutrally stable. As we 
go from a;* (2) towards the meeting point x m using type 2 
solutions, the integration is away from a node along the pri- 
mary direction and is unstable and extremely sensitive to the 
accuracy of the starting conditions. Based on intersections 
of the two curves in the a — v phase diagram, we identify 
matches for global similarity solutions. From these intersec- 
tion points, we determine the relevant £*(1) and x*(2) and 
compute solutions. Global similarity solutions obtained this 
way cross the magnetosonic critical line twice analytically. 

In Figure 1111 we present such a phase diagram of v 
versus a following the procedure described above. The first 
three intersection points for v > —0.5 in the phase di- 
agram lead to the first three global solutions that cross 
the magnetosonic critical line twice analytically. For the 
first solution, we have 0.10306 < x»(l) < 0.10307 and 
1.81092 < x»(2) < 1.81093. For the second solution, we have 
0.000128 < x t (l) < 0.000129 and 0.8729 < x«(2) < 0.8730. 
For the third solution, we have 1.2 x 10" 6 < < 
1.3 x 10 -6 and 1.165 < a;,(2) < 1.166. 

In Figure ED we display the first three global solutions 
that go across the magnetosonic critical line twice analyt- 
ically. Similar oscillation phenomena occur in magnetized 
cases as those in purely hydrodynamic cases llLou fc Shenl 
2004). In Figure ED we present the enlarged part in Figure 
1121 for small x with a logarithmic scale. 
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For A = 0, this equation reduces to the Lane-Emden equa- 
tion (see relevant discussions in the Appendix of Shu 1977) . 
These solutions represent self-similar inflows with a constant 
speed at large x. Besides matching with quasi-magnetostatic 
solutions, it is possible to match these solutions with self- 
similar outflows at smaller x with or without shocks (Tsai 
& Hsu 1995; Shu et al. 2002; Shen & Lou 2004; Bian & Lou 
2005). 

4-2.2 Self-Similar MHD Solutions crossing the 
Magnetosonic Critical Line Analytically 

To find the solutions crossing the magnetosonic critical line 
analytically, we follow the same procedure of Hunter (1977) 



4-2.3 Larson-Penston (LP) Type Solutions and Hunter 
Type Solutions 

iLarsonl jl969al'l and lPensto 3 |l969a) independently derived 
a kind of monotonic similarity solution that has the asymp- 
totic behaviour 1661 and 16711 with A = (i.e., without 
magnetic field) near the origin x — » + . In the presence 
of magnetic field, the MHD extension of LP solution has 
the asymptotic behaviour 1661 and 1671 with A > 0. The 
MHD LP-type solution can be obtained similarly by try- 
ing different values of D with a given A near the origin 
x — > + and by integrating nonlinear MHD ODEs 1241 and 
1251 rightward from x = towards the magnetosonic crit- 
ical line (x — v) + 2\x = (x — v) s until the critical condi- 
tions 1291 — 1311 are satisfied in such a way that the solution 
matches with one of the two eigensolutions across the mag- 
netosonic critical line. Thus the problem of finding global 
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Figure 11. The case of A = 0.1. The type 2-type 1 solutions 
matching phase diagram of v versus a. The meeting point is cho- 
sen at x m = 0.5. The type 2 eigenvalue is either positive in the v 
gradient (saddle point) or of larger absolute value for the v gradi- 
ent (nodal point). The type 1 is negative in the v gradient (saddle 
point) or of smaller absolute value for the v gradient (nodal point). 
The £C*(1) and x*(2) of the first solution of this type are in the 
range 0.10306 < x»(l) < 0.10307 and 1.81092 < x»(2) < 1.81093. 
For the second solution of this type, we have 0.000128 < x*(l) < 
0.000129 and 0.8729 < x*(2) < 0.8730. As for the third solu- 
tion of this type, we have 1.2 X 10~ 6 < x*(l) < 1.3 X 10~ 6 and 
1.165 < x»(2) < 1.166. 
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1st solution: x,(1)= 0.103 x,(2)=1.811 
2nd solution: X,(l)=l.28x10~* x,(2)=0.873 
3rd solution: x,(1)=1.2x10~ 6 x,(2)=1.165 
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Figure 12. The case of A = 0.1. The first three solutions passing 
across the magnetosonic critical line twice with increasing num- 
ber of stagnation points moving outward with sub-magnctosonic 
speeds. We have x»(l) « 0.10306 and x»(2) « 1.81092 for the 
first solution with one stagnation point; cc*(l) ~ 1.28 X 10 — 4 and 
x*(2) ~ 0.873 for the second solution with two stagnation points; 
and x»(l) »s 1.2 X 10~ 6 and x»(2) a 1.165 for the third solution 



with three stagnation points. The straight dashed line is x — 1 
and the dash-dotted curve is the magnetosonic critical line. 







LP-type solutions is basically an eigenvalue problem of de- 
termining discrete eigenvalues for the D parameter. 

For the hydrodynamic LP solution, the eigenvalue D 
is approximately 1.67. In our modified definition of type 1 
and type 2 solutions in the presence of magnetic field, the 
corresponding eigensolutions are related to type 1 solutions 



> 
I 




Log(x) 

Figure 13. The case of A = 0.1. The vertical coordinate is 
—v and the horizontal coordinate is log(x). By this enlargement, 
we see clearly that first the three solutions have one, two, and 
three stagnation points, respectively. The dash-dotted line is the 
magnetosonic critical line. 




Figure 14. The case of A = 0.1. The vertical coordinate is a 
and the horizontal coordinate is x. The first three solutions go 
across the magnetosonic critical line analytically. Solid, dashed 
and dash-dotted curves correspond to the first, second, and third 
solutions, respectively. 



[n.b., not type 2 solutions, see IShul il97?t): iHunterl Jl977t) : 
IWhitworth fe Summers! Jl985l) : iLou fe Shenl h004) }. 2 LP- 
type solution goes across the magnetosonic critical line an- 
alytically. From our numerical exploration for A = 0.1, 
the corresponding eigenvalue is D = 3.63. If we relax the 
constraints that solutions must be analytic (as defined by 
Hunter 1986) across the magnetosonic critical line, it is then 



IHunterl 11977 
I2004D . while the 



2 In previous studies JShu|_ 
IWhitworth fc Summersl Il985t lLou , 

LP solution is formally classified as type 2 solution (cf. our 
terminology) across the sonic critical line in the fourth quadrant, 
the corresponding eigensolution is actually along the secondary 
direction for x* > 2. Our type 1 and type 2 definitions would not 
suffer from the disadvantage of magnitude reversal in the sense of 
dv/dx. By our convention, type 1 is always along the secondary 
direction, and type 2 is always along the primary direction. 
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Figure 15. The case of A = 0.1. D = 2/(3 + 18A), 1.60, 2.40, 3.63 
running from lower to upper solution curves. LP-type solution 
passes across the magnetosonic critical line. The straight dashed 
line is x — v = 0. The slightly curved dash-dotted line is the 
magnetosonic critical line. 

possible to construct continuum bands of solutions cross- 
ing the magnetosonic critical line with weak discontinuities 
JWhitworth fc Summersil985l) . The case of D = 2/(3+ 18A) 
corresponds to the special exact global solution 1331 . 

Hunter type solutions refer to those similarity solutions 
crossing the magnetosonic critical line analytically and hav- 
ing asymptotic behaviour 16611 and 16711 as x approaches 0. 
These solutions pass through the magnetosonic critical line 
along the secondary direction, i.e., they pass the magne- 
tosonic critical line as type 1 solutions. We choose a meeting 
point at x m = 0.5. From around the origin x — > 0+ with D 
and A specified in equations 1661 and 167L we integrate the 
nonlinear MHD ODEs to the meeting point x m — 0.5. For 
every D, we obtain a {v, a} pair at x m in the phase diagram. 
By varying D parameter in sequence, we readily determine 
a series of {v, a} pairs. Plotting this series of {v, a}, we get 
a curve in the phase diagram. Meanwhile, we integrate the 
nonlinear MHD ODEs from different critical points x* (with 
x, > x m ) along the magnetosonic critical line backward to- 
wards the same meeting point x m to get another curve in 
the {a, v} phase diagram. Once the intersection points are 
determined in this phase diagram, we readily obtain the cor- 
responding Hunter type solutions. 



5 NOTES AND DISCUSSION 

In this paper, we investigated quasi-spherical self-similar 
flow solutions for magnetized self-gravitating isothermal flu- 
ids in the semi-complete solution space. The relevant MHD 
similarity solutions are obtained an d classified by compar- 
ing with results of previous analyses l|Larson| l969a|:|Penstoril 
1969d:IShJl977[lHunteJl977l:IWhitworth fc SummerJl985l 
Chiueh fe Choull994HLou fc SheiJ2004 IShen fc Loul2004 
The straight sonic critical line is now replaced by the curved 
magnetosonic critical line as expected physically. Global 
similarity solutions without crossing the magnetosonic criti- 
cal line also show ECCC and EECC characteristics, and are 
fairly easy to construct numerically in the isothermal ap- 
proximation. In a polytropic gas with radiative loss effects, 
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Figure 16. The case of A = 0.1. The meeting point is chosen at 
Xm = 0.5. The two parameters D and x t for the first solution of 
this type are in the ranges of 1.29 < In D < 1.30 and 2.164 < x* < 
2.165, respectively. For the second solution of this type, we have 
8.15 < In£> < 8.16 and 0.8704 < x* < 0.8705, respectively. And 
for the third solution of this type, we have 12.77 < InD < 12.78 
and 1.165 < x* < 1.166, respectively. 
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Figure 17. The case of A = 0.1 corresponding Fig. 1161 The first 
three Hunter type global solutions acrossing the magnetosonic 
critical line. The straight dashed line is x—v = 0. The dash-dotted 
curve is the magnetosonic critical line. For the first solution, we 
have lnZ5 = 1.29 and x* = 2.16; for the second solution, we have 
InD = 8.15 and x„ = 0.87; and for the third solution, we have 
InD = 12.7 and x» = 1.16. 

so me qualitatively similar feat ures were briefly mentioned 
in iBoilv fc Lvnden-Belll il995tl . It is also possible to con- 
struct similarity solutions for a magnetized polytropic gas 
with a quasi-spherical symmetry on larger scales (Wang & 
Lou 2005, in preparation). 

With magnetic field included, both magnetized expan- 
sion wave collapse solutions (mEWCS) and magnetized en- 
velope expansion core collapse (mEECC) solutions can be 
readily constructed. The mEWCS structure is composed 
of two main parts joined at the collapse wave front x» = 
(1 + 2A) 1/2 . The outer part with x > (1 + 2A) 1/2 is a nearly 
force-free magnetostatic envelope with the magnetic pres- 
sure and tension forces almost cancelling each other and 
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the inner part with x < (1 + 2A) 1 / 2 is an MHD core col- 
lapse characterized by diverging behaviours as x — » + . On 
the whole, such a mEWCS physically describes a situation 
with an MHD core collapse approaching a free-fall state in 
the very central region and with a magnetostatic exterior 
at rest for x ^ (1 + 2A) 1,/2 ; the quasi-spherical stagnation 
surface separating the interior and exterior expands at the 
magnetosonic speed (v\ + a 2 ) 1//2 . 

For more general mEECC similarity solutions with only 
one stagnation point x* along the x— axis, the interior col- 
lapse (a; < x*) and the exterior expansion [x ^ x*) concur in 
a self-similar manner. The outward expansion of the exterior 
envelope gives rise to a wind of constant speed at large radii. 
The quasi-spherical surface of separation where the gas re- 
mains at rest travels outward as time goes on, similar to the 
'expansion-wave front' of Shu (1977). From a more general 
perspective, the EWCS of Shu (1977) is a special case of our 
mEECC similarity solutions in the sense that an unmagne- 
tized outer envelope is static at large x instead of expanding 
or contracting at constant speeds. The travel speed of the 
stagnation front at x* could be either sub-magnetosonic for 
x* < (1 + 2A) 1 / 2 or super-magnetosonic for x* £ (1 + 2A) 1//2 . 

This explanation can be straightforward extended to 
those EECC solutions with more than one stagnation points, 
involving self-similar radial oscillations in sub-magnetosonic 
region. For such mEECC similarity solutions, there exist 
more than one quasi-spherical stagnation surfaces of sepa- 
ration that expand outward with time at different constant 
speeds. The separations between quasi-spherical surfaces in- 
crease with time in a self-similar manner. 

Such mEECC similarity solutions, which depict an out- 
going envelope and a collapsing core, are derived from non- 
linear MHD equations and may be applied to various as- 
trophysical processes with proper adaptions. For star form- 
ing regions in molecular clouds, for instance, the starless 
cloud system B335, we may take an infalling region of spa- 
tial scale ~ 1.5x 10 4 AU, a cloud mass of ~ 4 Mr , a magnetic 
field strength of ~ 134/iG dWolf et al]l2003ft and a thermal 
sound speed a ~ 0.23 km s _1 . In this case, the parameter A 
is approximately 0.01. By applying the EECC solution that 
crosses the magnetosonic critical line twice analytically with 
one stagnation point to model B335 system, we have the fol- 
lowing estimates. As cloud materials collapse from a radial 
distance 7700 AU to 2200 AU in our frozen-in scenario, the 
magnetic field strength at 2200 AU would become ~ 300/iG 
and the inflow speed at 2200AU is ~ 0.11 km s" 1 . This in- 
ward increase of 'primordial magnetic field' may serve as an 
important 'seed magnetic field' for a further enhancement 
through MHD dynamo processes within the central core of 
a protostar; the eventual appearance of generic bipolar out- 
flows which are most likely magnetized may result from ro- 
tation and symmetry breaking of gross magnetic field config- 
urations. For central collapse regions very much closer to the 
protostar, our self-similar solutions are not expected to be 
applicable as transients and feedbacks lead to complicated 
situations. 

Such mEECC similarity solutions with embedded mag- 
netic fields may also be applied to the asymptotic giant 
branch (AGB) phase or post-AGB phase in the late evolu- 
tion stage of a star before the gradual emergence evolution 
of a planetary nebula system with a central white dwarf of 
a high temperature. The timescale of this evolution 'gap' is 



estimated to be ~ 10 3 yrs. By the conventional wisdom, a 
main sequence star with a progenitor of less than several so- 
lar masses swells tremendously in size and sustains a massive 
wind with an asymptotic speed of ~ 10—20 km s - ; the mass 
loss rates fall within the range of ~ 10~ 8 to ~ 1O~ 4 M0 yr -1 . 
With an insufficient nuclear fuel supply from a certain epoch 
on, the central region starts to collapse while the outer enve- 
lope continues to expand; it is highly likely that the central 
collapse is accompanied by an outward energetic surge to 
chase the slowly expanding envelope. Simultaneously, the 
outer expansion removes stellar envelope mass, while the 
central infall and collapse produce a proto white dwarf at the 
center. Sufficiently far away from the initial and boundary 
conditions, the system may gradually evolve into a dynamic 
phase representable by an EECC similarity solution during a 
timescale of a few hundred to several thousand years. Within 
this rough range of timescale, the accumulation of central 
core mass should not exceed the Chandrasekhar mass limit 
of IAMq for a remnant white dwarf to survive; otherwise, 
a central explosion (e.g., novae or even a supernova) may 
occur within a planetary nebula. As a result of the energetic 
surge, the outer envelope expansion during the EECC phase 
may be faster than the pre-existing massive slow wind. This 
process will generally give rise to the formation of outward 
moving shock when a faster EECC outflow catches up with a 
slower pre-existing wind. Theoretically, it is possible to con- 
struct EECC similarity shock solutions (Shen & Lou 2004; 
Bian & Lou 2005) to generalize the earlier models of Tsai 
& Hsu (1995) and Shu et al. (2002). We therefore hypothe- 
size that dynamical evolution of an mEECC phase of around 
or less than a few thousand years may be the missing link- 
age between the AGB or post-AGB phase and the gradual 
appearance of a planetary nebula. Depending on physical 
parameters of the progenitor star, it might well happen that 
the core collapse and subsequent central infalls during an 
mEECC phase lead to an eventual core mass exceeding the 
Chandrasekhar mass limit of 1.4M© and thus induce an ex- 
plosion with an intensity determined by the actual rate of 
central mass accretion. 

For gross estimates, observations indicate a magnetic 
field strength of ~ 1G near the stellar surface at r ~ 1AU 
during the AGB phase, decreasing to ~ 10~ 3 G at r ~ 
10 3 AU. In the self-similar expansion regime of our model 
analysis, we have the transverse magnetic field £>x oc r _1 
(e.g., Lou 1993, 1994), consistent with the measured mag- 
netic field variations observed for AGB stars. Much stronger 
surface magnetic fields of a white dwarf might be generated 
and sustained by intrinsic MHD dynamo processes caused 
by convective differential rotation within the interior of an 
AGB star (e.g., Blackman et al. 2001). 

The transverse magnetic field strength of a mEECC so- 
lution in the innermost free-fall collapse region behaves as 
B±_ oc r -1 / 2 . For a magnetic field strength of 1G at r ~ 1AU 
and the radius of a proto-white dwarf of ~ 6000 km, the or- 
der of magnitude in spatial scale of the system decreases 
roughly by a factor of ~ 10 6 and accordingly the magnetic 
field strength increases roughly by a factor of ~ 10 3 . As 
a result, the magnetic field estimated at the surface of a 
proto-white dwarf is roug hly of the order of ~ 10 3 G. 

For the Crab Nebula as an example of supernova rem- 
nant, we simply model the magnetic field as B± oc r _1 dur- 
ing the envelope expansion phase, ignoring complex interac- 
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tions between the pulsar wind and the inner nebula involved 
in the expansion of the nebula. For a neutron star with a ra- 
dius of ~ 10 km and an initial dipolar magnetic field strength 
of ~ 10 10 G, the magnetic field would decrease to ~ 10 -3 G 
at about several parsecs. These estimates roughly agree with 
the envelope expansion portion of our simple mEECC solu- 
tions. 

For astrophysical systems of much larger scales such 
as clusters of galaxies (e.g., Sarazin 1988; Fabian 1994), 
the similarity solutions may be valuable for modeling 
a certain phase of their dynamical evolution involv- 
ing large-scale random magnetic fields. Specifically, if 
their evo lution involves a self-similar phase of central 
collapse <Gunn fe Gotd Il972t iFillmore fe Goldreichl Il984t 
iBertschingeJ Il985t iNavarro et alj Il997ft . then isothermal 
similarity solutions 16311 , 16411 , 1661 and I67H seem to suggest 
two possible classes of galaxy clusters that release X-ray 
emissions through hot electrons virialized in gravitational 
potential wells. The bifurcation of extremely high X-ray core 
luminosities and normal X-ray core luminosities may corre- 
spond to steeper and shallower gravitational potential wells, 
respectively. Typically, a galaxy cluster is characterized by 
a size of order of a few Mpc, a mass of order of 10 13 M Q , and 
a magnetic field strength of a few /uG. 

For a galaxy cluster resulting from a collapse of a grav- 
itationally bound system on much larger scales, a primodial 
weak magnetic field frozen into the gas can be enhanced by 
gradual contractions. In the central core collapse region, a 
transverse magnetic field may scale as B oc r _1//2 . In this 
scenario, specific values of primordial magnetic fields would 
be determined by the spatial scale when the system starts a 
self-similar collapse. 

While our mEECC and mECCC models are highly sim- 
plified and idealized, it would presumably grasp a few fun- 
damental aspects of certain physical processes and provide 
valuable physical insights into the magnetic field evolution 
in various astrophysical systems. 
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APPENDIX A: ANALYSIS OF SADDLE 
POINTS AND NODAL POINTS 

This appendix A is to show explicitly which portion of the 
magnetosonic critical line corresponds to saddle points and 
which portion of the magnetosonic critical line corresponds 
to nodal points (e.g., Jordan & Smith 1977). We examine 
equation I35|l in order to extract information about the signs 
of the roots of equation l|35|l along the magnetosonic critical 
line. In the hydrodynamical case, the roots of equation 1351 1 
can be expressed in fairly simple analytical forms. We know 
that a point within the open interval < x < 1 along the 
magnetosonic critical line is a saddle point, while a point 



within x > 1 along the magnetosonic critical line is a nodal 
point. 

In the presence of magnetic field, saddle and nodal 
points along the magnetosonic critical line are not so 
obvious as for the cases without magnetic field. Here, 
we perform such an analysis . Following the procedure of 
IWhitworth fc Summers! lll985l) , we cast the self-similar non- 



linear MHD ODEs in the forms of 

dv/dx = Y(x, v, a)/X (x, v, a) , (Al) 

da/dx = Z(x,v, a) /X(x,v, a) , (A2) 
where functionals X, Y, and Z are defined by 

X(x,v,a) = (x-v) 2 - (l + \ax 2 ) , (A3) 

Y(x, v, a) = a(x — v) 2 — 2(x — v)/x , (A4) 

Z(x,v,a) = a 2 (x — v) + 2\xa 2 — 2a(x — v) 2 /x . (A5) 



We figure out local topological properties of the solutions by 
examining the eigenvalues of the following matrix 

/ dX/dx dX/dv dX/da \ 
M= dY/dx dY/dv dY/da 

V dZ/dx dZ/dv dZ/da J 

with the partial derivatives of X, Y and Z evaluated at a 
given point along the magnetosonic critical line. Detailed 
calculations show that one of the three eigenvalues is zero 
and the other two eigenvalues satisfy the following quadratic 
equation 

az 2 + bz + c = , (A6) 

where a = 1 , b= -2£ , c = 2(C - x)(l - 3( 2 )/{{x 2 ) and C 
r r with £ > and x > under consideration. If the two 
roots of equation l|A6|l are of the same sign corresponding 
to the magnetosonic critical point being a nodal point, the 
requirement is simply c/a > 0, or equivalently ac > 0; this 
translates to (1 — 3£ 2 )(£ — x) > . For x < \/3/3 , we then 
have 

< -v < V3/3 - x . 

From Figure ijAlfl of the magnetosonic critical line, we see 
that this inequality cannot be realized. On the other hand 
for x > \/3/3 , it follows that 

73/3 - x < < . 

Points along the magnetosonic critical line in the range of 
x > %/3/3 which also satisfy the above inequality simulta- 
neously must be in the range of (1 + 2A) 1 / 2 < x < +oo . 

For the case of two real roots, the discriminant must be 
non-negative. The proof is as follows. In the open interval 
of (1 + 2A) 1 / 2 < x < +oo, we have the inequality f > (1 + 
2A) 1 / 2 > 1 and the discriminant is simply 

& 2 -4ac = 4C 2 -8^%^ 

_ 4x 2 ( 3 -8((-x)(l-3( 2 ) 

~ xH ' ( ' 

For a positive denominator, we only need to examine the 
sign of the numerator 

4z 2 C 3 -8(C-:r)(l-3C 2 ) 

= 4CV + (8 - 24< 2 )x - (8 - 24< 2 )< . 
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This quadratic expression in terms of x approaches the min- 
imum value of 

4(-l + 6C 2 -llC 4 + 6C 6 )/C 3 

as x takes on the value of x — (3£ 2 — 1)/C 3 > where the 
numerator can be further factorized into 

4(3C 2 - 1)(2C 2 - 1)(C + 1)(C ~ 1) ■ 

For £ > 1, this minimum remains always positive and thus 
the discriminant of the critical quadratic equation remains 
always positive. We thus come to the conclusion that when 
the abscissa x of a magnetosonic critical point falls within 
the open interval (1 + 2A) 1 / 2 < x < +oo, it is a nodal point. 

In the case of two real roots with opposite signs (i.e., 
a positive b 2 — 4ac with ac < 0) for a magnetosonic critical 
point being a saddle point, we have (1 — 3£ 2 )(£ — x) < . 
When x < s/E/3, we have either 



-x < -v < 



x < 



As seen from Figure ijAl^ . the latter is actually realizable. 
On the other hand when x > \/3/3, we have either 



< -v< V3/3 ■ 



< -v 



Again from Figure ijAlfl . the latter can be fulfilled along the 
magnetosonic critical line. Also by Figure (lAlfl . it is clear 
that the magnetosonic critical line always lies on the upper- 
right side of the straight line — v = \/3/3~ x, indicating — v > 
\/3/3 — x along the magnetosonic critical line. Combining 
these two analyzed cases for saddle type critical points, we 
come to the conclusion that for < x < (1 + 2A) 1//2 along 
the magnetosonic critical line, the critical point is of saddle 
type. 

We now summarize the key results of this Appendix A 
below. Along the magnetosonic critical line, for 

0<a;< (1 + 2A) 1/2 , 

the critical point is a saddle point, while for 

(1 + 2A) 1/2 < x < +oo , 

the critical point is a nodal point. 



APPENDIX B: DERIVATION OF THE 
QUADRATIC EQUATION FOR V' = Z 

Equation \Y2\t can be rearranged into the form of 

v = B/A , (Bl) 

where 

A = (x~v) 2 - (1 + Xax 2 ) , 
B ee (x - v)[a(x - v) - 2/x] . 

Along the magnetosonic critical line, both the denominator 
A and the numerator B of the above equation vanish. We 
use the L'Hospital rule to determine the derivatives at points 
along the magnetosonic critical line. 

The derivative of A with respect x is given by 

A' = 2(x - w)(l - v') - 2\ax - \x 2 a . 

From equation 12 1 1 . it follows that 
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Figure Al. Magnetic field parameter A = 1, the vertical coor- 
dinate is f = x — v, and the horizontal coordinate is x = r/(at). 
The S shaped curve is the magnetosonic critical line expressed in 
terms of variable f . Combining this figure with our analysis, we 
can readily distinguish nodal and saddle points along the magne- 
tosonic critical line. 



-(x — v) ■ 



(v — x) 



Using this expression in the first derivative of A, we know 
that A is a linear expression in terms of v = z, namely 

A = f(z) , 
where 



/(*) 



-2(a 



z + 2(x- 



(v — x) 

Meanwhile, we can express the first derivative of B as 

B = g(z) , 

where 

\_ 2«[1 — x(x — v)a] x[~ 2 + x(x — v)a] 
9(Z) = ~ 2 ~ 2 Z ■ 

Along the magnetosonic critical line, we have the critical 
condition that x(x — v)a — 2. With this relation, f(z) and 
g(z) can be reduced to the forms of 



/(*) 



2Aa; 
(x — v) 2 



- 2(x - v) 



z + 2(x- 



g(z) = 2v/x . 

By the L'Hospital rule, we know that z — g(z)/f(z) , or 
equivalently f(z)z — g(z) — , leading to quadratic equation 
ll?T51 l. 

Note that when A = and thus x — v = 1, we recover 
the results of Shu (1977), namely 

- 2z 2 + 2z - 2{x - 1) /x 2 = . (B2) 

The two roots of this equation are simply 

zi = l/s* and Z2 — 1 — 1/x, , 

respectively, where k* is a point along the sonic critical line. 
For higher-order expansion terms, see more details of Ap- 
pendix of Shu (1977) and of Lou & Shen (2004). 
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APPENDIX C: MAGNETIC PRESSURE AND 
MAGNETIC TENSION FORCES 

In spherical polar coordinates (r, 8, <j>), the three r— , 8— and 
(^—components of the magnetic tension force are: 



(B VB) 



B dBr B e dB r 
r dr r 88 



+ 



Ba, 8B r B 2 e+Bl 



(B ■ VB)s = Br^rr^- + 

or r 



r sm a o 
Be 8B e 



B<f, dBe Be B r 



cot 8Bl 



(B ■ VB) = B, 



r sin 8 d(f> 
SBj, , B e dB^ 



+ 



dr r d8 
B$ dBtf, B^B r 



+ 



cot 8B<j>Be 



(CI) 



(C2) 



(C3) 



Meanwhile, the radial component of the magnetic pressure 
gradient is 



^V(B 2 ) 



„ dB r 1 8 ,„j „o. 

B ^ + 2d-r^ + B ^ 

-Br d -^ + Be d -^ + B, d -§± 
dr dr dr 



(C4) 



the 8 component of the magnetic pressure gradient is 



^V(B 2 ) 



B r dB r Be dBe B$ dB^ 



(C5) 



/ kj \J I \J \J I 1/1/ 

and the 4> component of the magnetic pressure gradient is 



-V(B 2 



2r sin ( 



-{Bl + B 2 e + Bl) 



B r dB r 



+ 



Be dB e 



+ ■ 



B s dB. 







8 d<t> ' y ' 

Strictly speaking, when the magnetic field is taken into ac- 
count, a global spherical symmetry cannot be sustained. 
Checking through equations IC1I . 1C2L 1C3II . and equations 
1C4I I , 1C5I I , lC6t , we see that second and third terms in equa- 
tion iCll are ignored. The first term in (ICH and the first 
term in (IC4t cancel each other. In the radial momentum 
equation, the fourth term in equation ICH and the second 
and third terms in equation 1C41 are present and kept. 

As we take the magnetofluid system to be quasi- 
spherically symmetric on larger scales, the terms involving 
d/d8 and d/d(/> are neglected. Thus in the 8 and <j> direc- 
tions, the first, fourth and fifth terms in equations <C2t and 
dCJ3l > would cause the overall magnetofluid to be nonspher- 
ically symmetric. In order to work in the framework of a 
quasi-spherical symmetry, we would drop these these terms. 
Nevertheless, we should keep in mind that these terms might 
give rise to nonspherical shapes of various nebulae in con- 
texts of astrophysics. The quasi-elliptical morphology of the 
Crab Nebula is a concrete example. 



